#ifndef SPLINE_H
#define SPLINE_H

#include <cblas.h> // 对于BLAS
#include <lapacke.h> // 对于LAPACK
#include <stdexcept>
#include <cmath> 
#include <vector>
#include <iostream>

class Spline {
protected:
    const std::vector<double> node;
public:
    Spline(const std::vector<double>& node) : node(node) {}

    virtual double SplineReturn(double x) const = 0; 

    virtual ~Spline() = default;
};

class PPFormSpline : public Spline {
private:
    std::vector<double> first_derivatives;  // 存储每个节点的一阶导数
    std::vector<double> y;  // 存储样条节点的函数值
    std::vector<double> dividedDifferences;  // 用于存储差商的中间计算结果

    // 成员函数：将数组 b 复制到 first_derivatives
    void copyToFirstDerivatives(double b[], int n) {
        first_derivatives = std::vector<double>(b, b + n);
    }
    
    // 计算三对角方程的解
    static std::vector<double> solveTridiagonal(
        const std::vector<double>& lower,  // 下三角元素
        const std::vector<double>& diag,   // 对角线元素
        const std::vector<double>& upper,  // 上三角元素
        const std::vector<double>& right   // 右边的常数项
    ) {
        int n = diag.size();
        std::vector<double> c(n), l(n), u(n);
        
        // 初始化第一个方程
        l[0] = diag[0];
        u[0] = upper[0] / l[0];
        c[0] = right[0] / l[0];
        
        // 前向消元
        for (int i = 1; i < n; ++i) {
            l[i] = diag[i] - lower[i-1] * u[i-1];
            u[i] = upper[i] / l[i];
            c[i] = (right[i] - lower[i-1] * c[i-1]) / l[i];
        }
        
        // 回代
        std::vector<double> result(n);
        result[n-1] = c[n-1];
        for (int i = n-2; i >= 0; --i) {
            result[i] = c[i] - u[i] * result[i+1];
        }
        
        return result;
    }

    // 自然边界条件
    void computeNaturalBoundary(const std::vector<double>& x, const std::vector<double>& y) {
        int n = x.size();
        std::vector<double> lower(n-1), upper(n-1), diag(n, 2.0), right(n, 0.0);

        // 自然样条，边界条件：s''(a)=0,s''(b)=0
        lower[n-2] = 1.0;
        for (int i = 0; i < n-2; ++i) {
            lower[i] = (x[i+3]-x[i+2])/(x[i+3]-x[i+1]);
        }
        upper[0] = 1.0;
        for (int i = 1; i < n-1; ++i) {
            upper[i] = (x[i+1]-x[i])/(x[i+2]-x[i]);
        }
        right[0] = 3*dividedDifferences[0];  
        right[n-1] = 3*dividedDifferences[n-2]; 
        for (int i = 1; i < n-1; ++i) {
            right[i] = 3*upper[i]*dividedDifferences[i]+3*lower[i-1]*dividedDifferences[i-1];
        }

        // 求解方程得到一阶导数
        first_derivatives = solveTridiagonal(lower, diag, upper, right);
    }

    // 周期边界条件
    void computePeriodicBoundary(const std::vector<double>& x, const std::vector<double>& y) {
        int n = x.size();
        double A[n*n];  // 存储系数矩阵 A，大小为 n*n
        double b[n];   // 存储常数向量 b，大小为 n
        std::fill(A, A + n * n, 0.0);
        std::fill(b, b + n, 0.0);
        int ipiv[n];   // 存储主元素的排列
        int lda = n;    // A 的 leading dimension
        int ldb = 1;    // b 的 leading dimension
        int info;       // LAPACKE 返回的状态

        // 周期样条，边界条件：s(a)=s(b),s'(a)=s'(b),s''(a)=s''(b)
        A[0] = 2 / (x[1] - x[0]);
        A[1] = 1 / (x[1] - x[0]);
        A[n-2] = 1 / (x[n-1] - x[n-2]);
        A[n-1] = 2 / (x[n-1] - x[n-2]);

        for (int i = 1; i < n - 1; ++i) {
            A[i * n + (i - 1)] = (x[i + 1] - x[i]) / (x[i + 1] - x[i - 1]);
            A[i * n + i] = 2;
            A[i * n + (i + 1)] = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
        }

        // 处理周期边界条件
        A[n * n - n] = 1;
        A[n * n - 1] = -1;

        // 填充常数向量 b
        b[0] = 3 * dividedDifferences[0] / (x[1] - x[0]) + 3 * dividedDifferences[n - 2] / (x[n - 1] - x[n - 2]);
        for (int i = 1; i < n - 1; ++i) {
            b[i] = 3 * (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]) * dividedDifferences[i]
                + 3 * (x[i + 1] - x[i]) / (x[i + 1] - x[i - 1]) * dividedDifferences[i - 1];
        }   

        // 使用 LAPACKE 求解线性方程组 Ax = b
        info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, 1, A, lda, ipiv, b, ldb);

        // 打印解向量 b
        if (info == 0) {
            copyToFirstDerivatives(b,n);
        } else {
            std::cerr << "Error in solving the system, info = " << info << std::endl;
        }
    }

    // 完全边界条件
    void computeCompleteBoundary(const std::vector<double>& x, const std::vector<double>& y, const double fa, const double fb) {
        int n = x.size();
        std::vector<double> lower(n-3), upper(n-3), diag(n-2, 2.0), right(n-2, 0.0);

        // 完全样条，边界条件：s'(a)=f'(a),s'(b)=f'(b)
        for (int i = 0; i < n-3; ++i) {
            lower[i] = (x[i+3]-x[i+2])/(x[i+3]-x[i+1]);
        }
        for (int i = 0; i < n-3; ++i) {
            upper[i] = (x[i+1]-x[i])/(x[i+2]-x[i]);
        }
        right[0] = 3*upper[0]*dividedDifferences[1]+3*(x[2]-x[1])/(x[2]-x[0])*dividedDifferences[0]-(x[2]-x[1])/(x[2]-x[0])*fa;  
        right[n-3] = 3*(x[n-2]-x[n-3])/(x[n-1]-x[n-3])*dividedDifferences[n-2]+3*lower[n-4]*dividedDifferences[n-3]-(x[n-2]-x[n-3])/(x[n-1]-x[n-3])*fb; 
        for (int i = 1; i < n-3; ++i) {
            right[i] = 3*upper[i]*dividedDifferences[i+1]+3*lower[i-1]*dividedDifferences[i];
        }

        // 求解方程得到一阶导数
        first_derivatives = solveTridiagonal(lower, diag, upper, right);
        first_derivatives.insert(first_derivatives.begin(), fa);
        first_derivatives.push_back(fb);
    }

    // 计算一阶导数
    void computeFirstDerivatives(const std::vector<double>& x, const std::vector<double>& y, const std::string& boundary_type, const double fa, const double fb) {
        // 根据边界类型调用相应的函数
        if (boundary_type == "natural") {
            computeNaturalBoundary(x, y);
        } else if (boundary_type == "periodic") {
            computePeriodicBoundary(x, y);
        } else if (boundary_type == "complete") {
            computeCompleteBoundary(x, y, fa, fb);
        }
    }

    void computeFirstDerivatives(const std::vector<double>& x, const std::vector<double>& y, const std::string& boundary_type) {
        // 根据边界类型调用相应的函数
        if (boundary_type == "natural") {
            computeNaturalBoundary(x, y);
        } else if (boundary_type == "periodic") {
            computePeriodicBoundary(x, y);
        } 
    }

public:
    PPFormSpline(const std::vector<double>& x, const std::vector<double>& y, const std::string& boundary_type, const double fa, const double fb)
        : Spline(x), y(y) {
        int n = x.size();
        dividedDifferences.resize(n-1);
        
        // 计算差商
        for (int i = 0; i < n-1; ++i) {
            dividedDifferences[i] = (y[i+1] - y[i]) / (x[i+1] - x[i]);
        }

        // 计算一阶导数
        computeFirstDerivatives(x, y, boundary_type, fa, fb);
    }

    PPFormSpline(const std::vector<double>& x, const std::vector<double>& y, const std::string& boundary_type)
        : Spline(x), y(y) {
        int n = x.size();
        dividedDifferences.resize(n-1);
        
        // 计算差商
        for (int i = 0; i < n-1; ++i) {
            dividedDifferences[i] = (y[i+1] - y[i]) / (x[i+1] - x[i]);
        }

        // 计算一阶导数
        computeFirstDerivatives(x, y, boundary_type);
    }

     // 插值函数：根据一阶导数进行插值
    double SplineReturn(double x) const override {
        int n = node.size();
        int i = 0;
        
        // 查找区间所在的节点区间
        for (i = 0; i < n - 1; ++i) {
            if (x >= node[i] && x <= node[i+1]) break;
        }
        
        // 计算样条插值
        double h = node[i+1] - node[i];

        // 使用一阶导数进行插值
        double a = y[i];
        double b = first_derivatives[i];
        double c = (3*dividedDifferences[i]-2*first_derivatives[i]-first_derivatives[i+1])/(node[i+1]-node[i]);
        double d = (first_derivatives[i]+first_derivatives[i+1]-2*dividedDifferences[i])/std::pow(node[i+1]-node[i], 2);
        
        return a + b * (x - node[i]) + c * std::pow(x - node[i],2) + d * std::pow(x - node[i],3);
    }
};

class BSpline : public Spline {
private:
    int degree;
    std::vector<double> points;
    std::vector<double> y;
    std::vector<double> a;

    void copyToa(double b[], int n) {
        a = std::vector<double>(b, b + n);
    }

    // 自然边界条件
    void computeNaturalBoundary() {
        int m = points.size();
        int n = m-1-degree;
        double A[n*n];  // 存储系数矩阵 A，大小为 n*n
        double b[n];   // 存储常数向量 b，大小为 n
        std::fill(A, A + n * n, 0.0);
        std::fill(b, b + n, 0.0);
        int ipiv[n];   // 存储主元素的排列
        int lda = n;    // A 的 leading dimension
        int ldb = 1;    // b 的 leading dimension
        int info;       // LAPACKE 返回的状态

        // 自然样条，边界条件：s''(a)=0,s''(b)=0
        for (int i = 0; i < n; ++i) {
            A[i] = B_derivative(i,degree,points[degree],2);
        }

        for (int i = 1; i < n-1; ++i) {
            for (int j = 0; j < n; ++j)
            A[i*n+j] = B(j,degree,points[degree+i-1]);

        }

        for (int i = 0; i < n; ++i) {
            A[n*n-n+i] = B_derivative(i,degree,points[n],2);
        } 

        // 填充常数向量 b
        for (int i = 1; i < n-1; ++i) {
            b[i] = y[i-1];
        } 

        // 使用 LAPACKE 求解线性方程组 Ax = b
        info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, 1, A, lda, ipiv, b, ldb);

        // 打印解向量 b
        if (info == 0) {
            copyToa(b,n);
        } else {
            std::cerr << "Error in solving the system, info = " << info << std::endl;
        }
    }

    // 周期边界条件
    void computePeriodicBoundary() {
        int m = points.size();
        int n = m-1-degree;
        double A[n*n];  // 存储系数矩阵 A，大小为 n*n
        double b[n];   // 存储常数向量 b，大小为 n
        std::fill(A, A + n * n, 0.0);
        std::fill(b, b + n, 0.0);
        int ipiv[n];   // 存储主元素的排列
        int lda = n;    // A 的 leading dimension
        int ldb = 1;    // b 的 leading dimension
        int info;       // LAPACKE 返回的状态

        // 周期样条，边界条件：s(a)=s(b),s'(a)=s'(b),s''(a)=s''(b)
        for (int i = 0; i < n; ++i) {
            A[i] = B_derivative(i,degree,points[degree],2)-B_derivative(i,degree,points[n],2);
        }

        for (int i = 1; i < n-1; ++i) {
            for (int j = 0; j < n; ++j)
            A[i*n+j] = B(j,degree,points[degree+i-1]);

        }

        for (int i = 0; i < n; ++i) {
            A[n*n-n+i] = B_derivative(i,degree,points[degree],1)-B_derivative(i,degree,points[n],1);
        } 

        // 填充常数向量 b
        for (int i = 1; i < n-1; ++i) {
            b[i] = y[i-1];
        } 

        /*// 打印矩阵 A
        std::cout << "Matrix A:" << std::endl;
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                std::cout << A[i * n + j] << " "; // 打印矩阵元素
            }
            std::cout << std::endl;
        }

        // 打印解向量 b
        std::cout << "Solution to the system (first derivatives):" << std::endl;
        for (int i = 0; i < n; ++i) {
            std::cout << b[i] << std::endl; // 打印向量 b
        }*/

        // 使用 LAPACKE 求解线性方程组 Ax = b
        info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, 1, A, lda, ipiv, b, ldb);

        // 打印解向量 b
        if (info == 0) {
            copyToa(b,n);
        } else {
            std::cerr << "Error in solving the system, info = " << info << std::endl;
        }
    }

    // 完全边界条件
    void computeCompleteBoundary(const double fa, const double fb) {
        int m = points.size();
        int n = m-1-degree;
        double A[n*n];  // 存储系数矩阵 A，大小为 n*n
        double b[n];   // 存储常数向量 b，大小为 n
        std::fill(A, A + n * n, 0.0);
        std::fill(b, b + n, 0.0);
        int ipiv[n];   // 存储主元素的排列
        int lda = n;    // A 的 leading dimension
        int ldb = 1;    // b 的 leading dimension
        int info;       // LAPACKE 返回的状态

        // 完全样条，边界条件：s'(a)=f'(a),s'(b)=f'(b)
        for (int i = 0; i < n; ++i) {
            A[i] = B_derivative(i,degree,points[degree],1);
        }

        for (int i = 1; i < n-1; ++i) {
            for (int j = 0; j < n; ++j)
            A[i*n+j] = B(j,degree,points[degree+i-1]);

        }

        for (int i = 0; i < n; ++i) {
            A[n*n-n+i] = B_derivative(i,degree,points[n],1);
        } 

        // 填充常数向量 b
        b[0]=fa;
        for (int i = 1; i < n-1; ++i) {
            b[i] = y[i-1];
        } 
        b[n-1]=fb;

        /*// 打印矩阵 A
        std::cout << "Matrix A:" << std::endl;
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                std::cout << A[i * n + j] << " "; // 打印矩阵元素
            }
            std::cout << std::endl;
        }

        // 打印解向量 b
        std::cout << "Solution to the system (first derivatives):" << std::endl;
        for (int i = 0; i < n; ++i) {
            std::cout << b[i] << std::endl; // 打印向量 b
        }*/

        // 使用 LAPACKE 求解线性方程组 Ax = b
        info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, 1, A, lda, ipiv, b, ldb);

        // 打印解向量 b
        if (info == 0) {
            copyToa(b,n);
        } else {
            std::cerr << "Error in solving the system, info = " << info << std::endl;
        }
    }

    void compute_a(const std::string& boundary_type, const double fa, const double fb) {
        // 根据边界类型调用相应的函数
        if (boundary_type == "natural") {
            computeNaturalBoundary();
        } else if (boundary_type == "periodic") {
            computePeriodicBoundary();
        } else if (boundary_type == "complete") {
            computeCompleteBoundary(fa, fb);
        }
    }

    void compute_a(const std::string& boundary_type) {
        // 根据边界类型调用相应的函数
        if (boundary_type == "natural") {
            computeNaturalBoundary();
        } else if (boundary_type == "periodic") {
            computePeriodicBoundary();
        } 
    }

    void compute_a() {
        int m = points.size();
        int n1 = m-1-degree; // A的列数
        int n2 = m-2*degree; // A的行数
        double A[n2*n1];  // 存储系数矩阵 A，大小为 n*n
        double b[n2];   // 存储常数向量 b，大小为 n
        std::fill(A, A + n2*n1, 0.0);
        std::fill(b, b + n2, 0.0);
        int lda = n1;    // A 的 leading dimension
        int ldb = 1;    // b 的 leading dimension
        int info;       // LAPACKE 返回的状态

        for (int i = 0; i < n2; ++i) {
            for (int j = 0; j < n1; ++j)
            A[i*n1+j] = B(j,degree,points[degree+i]);
        }

        // 填充常数向量 b
        for (int i = 0; i < n2; ++i) {
            b[i] = y[i];
        } 

        // 使用 LAPACKstd::coutE 求解欠定线性方程组 Ax = b
        info = LAPACKE_dgels(LAPACK_ROW_MAJOR, 'N', n2, n1, 1, A, lda, b, ldb);

        // 打印解向量 b
        if (info == 0) {
            copyToa(b,n1);
        } else {
            std::cerr << "Error in solving the system, info = " << info << std::endl;
        }
    }

    // 计算B样条基函数
    double B(int i, int degree, double t) const {
        if (degree == 0) {
            if (points[i] <= t && t < points[i+1]) return 1.0;
            return 0.0;
        }
        double left = (t - points[i]) / (points[i+degree] - points[i]);
        double right = (points[i+degree+1] - t) / (points[i+degree+1] - points[i+1]);
        return left * B(i, degree-1, t) + right * B(i+1, degree-1, t);
    }

    // 计算B样条m阶导数
    double B_derivative(int i, int degree, double t, int m) const {
        double left = degree / (points[i+degree] - points[i]);
        double right = degree / (points[i+degree+1] - points[i+1]);
        if (m == 0) {
            return B(i, degree, t);
        }
        
        return left * B_derivative(i, degree-1, t, m-1) - right * B_derivative(i+1, degree-1, t, m-1);
    }

public:
    BSpline(const std::vector<double>& x, const std::vector<double>& y, int degree, const std::string& boundary_type, const double fa, const double fb)
        : Spline(x), y(y), degree(degree), points(x) {
            int N = node.size();
            double L = node[N-1]-node[0];
            for (int i = 0; i < degree; ++i){
                points.insert(points.begin(), node[0]-(i+1)*L/(N-1));
                points.push_back(node[N-1]+(i+1)*L/(N-1));
            }
            if (degree==3){
                compute_a(boundary_type, fa, fb);
            }
        }

    BSpline(const std::vector<double>& x, const std::vector<double>& y, int degree, const std::string& boundary_type)
        : Spline(x), y(y), degree(degree), points(x) {
            int N = node.size();
            double L = node[N-1]-node[0];
            for (int i = 0; i < degree; ++i){
                points.insert(points.begin(), node[0]-(i+1)*L/(N-1));
                points.push_back(node[N-1]+(i+1)*L/(N-1));
            }
            if (degree==3){
                compute_a(boundary_type);
            }
        }

    BSpline(const std::vector<double>& x, const std::vector<double>& y, int degree)
        : Spline(x), y(y), degree(degree), points(x){
            int N = node.size();
            double L = node[N-1]-node[0];
            for (int i = 0; i < degree; ++i){
                points.insert(points.begin(), node[0]-(i+1)*L/(N-1));
                points.push_back(node[N-1]+(i+1)*L/(N-1));
            }
            compute_a();
        }

    BSpline(const std::vector<double>& x, const std::vector<double>& y, int degree, const std::vector<double>& a)
        : Spline(x), y(y), degree(degree), points(x),a(a) {
            int N = node.size();
            double L = node[N-1]-node[0];
            for (int i = 0; i < degree; ++i){
                points.insert(points.begin(), node[0]-(i+1)*L/(N-1));
                points.push_back(node[N-1]+(i+1)*L/(N-1));
            }
        }

    double SplineReturn(double x) const override {
        int m = points.size();
        int n = m-1-degree;
        double result = 0.0;
        
        // 计算基函数B的值，并加权
        for (int i = 0; i < n; ++i) {
            result += a[i] * B(i, degree, x);
        }
        return result;
    }
};

#endif